R Notebook

data prep

library(Seurat)
library(Matrix)

combined_cell@meta.data$cell_type <- Idents(combined_cell)
combined_cell@meta.data$organ_celltype <- Idents(combined_celltype)
combined_cell@meta.data$organ <- Idents(combined_organ)

# Load Seurat object
options(future.globals.maxSize = 2000 * 1024^2)
Idents(combined_cell) <- combined_cell@meta.data$organ
all.markers <- FindAllMarkers(object = combined_cell,verbose = FALSE, logfc.threshold = 0.2)
write.csv(all.markers, "deg_across_organs.csv")

brain.pericyte.markers <- FindMarkers(object = combined_cell, ident.1 = "Brain_mural", ident.2 = c("Bladder_mural", "Heart_mural", "Colon_mural", "Muscle_mural"), group.by = "organ_celltype", verbose = FALSE, logfc.threshold = 0.2)

For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
(default method for FindMarkers) please install the presto package
--------------------------------------------
install.packages('devtools')
devtools::install_github('immunogenomics/presto')
--------------------------------------------
After installation of presto, Seurat will automatically use the more 
efficient implementation (no further action necessary).
This message will be shown once per session

write.csv(brain.pericyte.markers, "deg_pericyte.csv")

write.csv(as.matrix(combined_cell[,1:2000]@assays$RNA@counts), "expression_matrix_test.csv")
write.csv(combined_cell[,1:2000]@meta.data, "metadata_test.csv")
write.csv(combined_cell[,1:2000]@assays$RNA@features, "gene_features_test.csv")

Error in is.data.frame(x) : 
  no slot of name "features" for this object of class "Assay"

working dir

setwd("~/Library/CloudStorage/GoogleDrive-schen601@usc.edu/My Drive/transportor/transporter_bert")

deg filtering

# break marker table
library(dplyr)

############### split degs: all.markers ###############
# create an output directory to store the files
output_dir <- file.path(dirname(getwd()), "deg")
dir.create(output_dir, showWarnings = FALSE)

# Split data by the "organ" column and save each subset
split_data <- split(all.markers, all.markers$cluster)

# Save each subset as a separate CSV file
for (organ_name in names(split_data)) {
  file_name <- paste0(output_dir, "/", organ_name, ".csv")
  write.csv(split_data[[organ_name]], file_name, row.names = FALSE)
}

# - ignore
slc_map <- read_csv("group-752(slc_family_HGNC).csv")
gene_symbols <- rownames(combined_cell[["RNA"]])
gene_symbols <- unique(gene_symbols)

library(tidyr)
required_cols <- c("Approved symbol", "Previous symbols", "Aliases", "Alias symbols", "Approved.symbol", "Previous.symbols")

# Convert the relevant columns to character
slc_map <- slc_map %>%
  mutate(across(all_of(required_cols), as.character))

# genes from SLC family
gene_symbols_slc <- slc_map %>%
  select(all_of(required_cols)) %>%
  pivot_longer(cols = everything(), values_to = "gene") %>%
  separate_rows(gene, sep = ",\\s*") %>%   # Split by comma and remove spaces
  distinct() %>%  # Remove duplicate entries
  na.omit()       # Remove NA values
gene_symbols_slc <- unique(gene_symbols_slc$gene)


# shared slc
slc_targets <- intersect(tolower(gene_symbols), tolower(gene_symbols_slc))

# filter slc degs
filtered.markers <- all.markers %>%
  mutate(gene = tolower(gene)) %>% 
  filter(gene %in% slc_targets)  # Keep only rows where gene is in slc_targets

############### split degs: filetred.markers ###############
# create an output directory to store the files
output_dir <- file.path(dirname(getwd()), "deg_SLC")
dir.create(output_dir, showWarnings = FALSE)

# Split data by the "organ" column and save each subset
split_data <- split(filtered.markers, filtered.markers$cluster)

# Save each subset as a separate CSV file
for (organ_name in names(split_data)) {
  file_name <- paste0(output_dir, "/", organ_name, ".csv")
  write.csv(split_data[[organ_name]], file_name, row.names = TRUE)
}

dot plot

# inputs HGNC dataset -> save into c(df)
hgnc_dir <- "/Users/chen/Library/CloudStorage/GoogleDrive-schen601@usc.edu/My Drive/transportor/hgnc/"
hgnc_filename_list <- list.files(path = hgnc_dir, full.names = TRUE)

# check file type and read accordingly
read_hgnc_smart <- function(file) {
  first_line <- readLines(file, n = 1, warn = FALSE)
  # Check if the first line contains "sep=,"
  if (grepl("^sep=,", first_line)) {
    return(read.csv(file, skip = 1))  # Read as CSV, skipping first line
  } else {
    return(read_delim(file, delim = "\t", escape_double = FALSE, trim_ws = TRUE))
  }
}

hgnc_df_list <- lapply(hgnc_filename_list, read_file_smart)
names(hgnc_df_list) <- basename(hgnc_filename_list)

generate_dotplot <- function(combined_cell, sorted_feature_genes, file_name) {
  dot_plot <- DotPlot(combined_cell, features = rev(sorted_feature_genes), group.by = "organ_sorted") +
  labs(x = "Organ", y = "Gene") +
  scale_x_discrete(position = "bottom") +  # Move organ labels to the top
  scale_y_discrete(position = "right") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        axis.text.x.top = element_text(angle = 45, hjust = 0, vjust = 0),
        axis.ticks.x.top = element_line(),
        axis.title.x = element_blank()) +
  ggtitle(gsub("\\.csv$", "", file_name)) +
  coord_flip() 
  print(dot_plot)
  
  # save into phf
  plot_height <- max(2, (1.2 + (length(sorted_feature_genes) * 0.3)))
  ggsave(paste0("dot_plot_", gsub("\\.csv$", "", file_name), ".pdf"), plot = dot_plot, width = 7, height = plot_height,limitsize = FALSE)
}

sort_gene_by_ubiquity <- function(combined_cell, feature_genes, file_name) {
  # Extract expression data for the selected genes
  expr_matrix <- as.data.frame(GetAssayData(combined_cell, assay = "RNA", slot = "data")[feature_genes, ])
  #if (file_name == "Sodium_channels.csv"){
  if (length(feature_genes) == 1){
    expr_matrix <- as.data.frame(t(expr_matrix))
    row.names(expr_matrix) <- c(feature_genes)
  }

  # Convert matrix to data frame and add organ labels
  expr_df <- as.data.frame(expr_matrix) %>% rownames_to_column(var = "gene") %>% pivot_longer(cols = -gene, names_to = "cell", values_to = "expression")
  # Ensure metadata row names are accessible
  meta_data_df <- combined_cell@meta.data %>%
    rownames_to_column(var = "cell") %>%  # Convert row names to "cell" column
    select(cell, organ)  # Select relevant columns
  # Perform left join with metadata
  expr_df <- expr_df %>%
    left_join(meta_data_df, by = "cell")
  
  # Count the number of unique organs in which each gene is expressed
  gene_organ_count <- expr_df %>%
    group_by(gene) %>%
    filter(expression > 0) %>%  # Consider only expressed genes
    summarize(num_organs = n_distinct(organ)) %>%
    arrange(desc(num_organs))  # Sort by ubiquity
  
  # return genes
  return(gene_organ_count$gene)
}

# intersection with HGNC
#output_dir_hgnc_filtered_deg <- file.path(dirname(getwd()), "deg_brain_filtered_by_HGNC")
output_dir_hgnc_filtered_deg <- file.path(dirname(getwd()), "pericyte_DEG_filtered_by_HGNC")
dir.create(output_dir_hgnc_filtered_deg, showWarnings = FALSE)

# each transporter family
transporter_dictionary <- list()
filenames <- names(hgnc_df_list)
#filenames <- c("Glycine_receptors.csv")
for (file_name in filenames) {
  
  # =====================intersection: HGNC & brain========================
  # Get the current table
  trasporter_map_i <- hgnc_df_list[[file_name]]
  
  # Convert the interest columns to character
  exist_cols <- colnames(trasporter_map_i)[colnames(trasporter_map_i) %in% required_cols]
  trasporter_map_i <- trasporter_map_i %>%
    mutate(across(all_of(exist_cols), as.character))
  
  # genes from HGNC
  gene_symbols_hgnc_i <- trasporter_map_i %>%
    select(all_of(exist_cols)) %>%
    pivot_longer(cols = everything(), values_to = "gene") %>%
    separate_rows(gene, sep = ",\\s*") %>%   # Split by comma and remove spaces
    distinct() %>%  # Remove duplicate entries
    na.omit()       # Remove NA values
  gene_symbols_hgnc_i <- unique(gene_symbols_hgnc_i$gene)

  # shared slc
  
  #transporter_targets <- intersect(tolower(brain.markers.sandra$Gene), tolower(gene_symbols_hgnc_i))
  brain.pericyte.markers$Gene <- row.names(brain.pericyte.markers)
  transporter_targets <- intersect(tolower(brain.pericyte.markers$Gene), tolower(gene_symbols_hgnc_i))

  # Filter the table based on matched genes
  #filtered_table <- brain.markers.sandra %>%
  #  filter(tolower(Gene) %in% transporter_targets)
  filtered_table <- brain.pericyte.markers %>%
    filter(tolower(Gene) %in% transporter_targets)
  
  # Restore original case for gene names
  #filtered_table$gene <- brain.markers.sandra$Gene[brain.markers.sandra$gene %in% transporter_targets]
  
  print(gsub("\\.csv$", "", file_name))
  if (length(filtered_table$Gene) == 0) {
    next
  }
  cat(paste(c(length(filtered_table$Gene), "genes: ",paste0(filtered_table$Gene, " "))))

  # re-save the new DEGs (shared)
  output_file <- file.path(output_dir_hgnc_filtered_deg, paste0("filtered_brain_pc_DEGs_by_", file_name))
  write.csv(filtered_table, output_file, row.names = FALSE)
  
  # =====================prepare plot inputs: organs and genes========================
  if (length(filtered_table$Gene) < 25) {
    feature_genes <- filtered_table$Gene
  } else {
    feature_genes <- filtered_table$Gene[1:25]
  }
  
  # Brain sorted to top
  combined_cell@meta.data$organ_sorted <- factor(combined_cell@meta.data$organ, levels = c("Brain", setdiff(unique(combined_cell@meta.data$organ), "Brain")))

  # sort genes by ubiquiousness
  sorted_feature_genes <- sort_gene_by_ubiquity(combined_cell, feature_genes, file_name)
  transporter_dictionary[[gsub("\\.csv$", "", file_name)]] <- sorted_feature_genes
  
  # ================================dot plot===================================
  generate_dotplot(combined_cell, sorted_feature_genes, file_name)
  }

[1] "ABC_transporters"
2 genes:  Abcc9  Abca9 [1] "Amine_receptors"
[1] "ATPase"
7 genes:  Atp13a5  Atp1a2  Atp1b2  Atp2a3  Atp2c1  Atp10d  Atp7a [1] "Calcium_channels"
5 genes:  Itpr2  Cacna1c  Tpcn1  Itpr1  Cacna1h [1] "Cholinergic_receptors"
[1] "GABA_receptors"
[1] "Glutamate_receptors"
3 genes:  Grm3  Grm7  Glud1 [1] "Glycine_receptors"
1 genes:  Glrb [1] "GPCR"
16 genes:  P2ry14  Gper1  Grm3  S1pr3  Cd97  Pth1r  Tbxa2r  Grm7  Gprc5c  Fzd6  Adora2a  Gpr116  Gpr124  Ednra  Gpr4  Ackr3 [1] "Potassium_channels"
1 genes:  Kcnj8 [1] "SLC"
25 genes:  Slc22a8  Slc30a10  Slc19a1  Slc16a12  Slc12a2  Slc6a13  Slc5a5  Slc38a11  Slc25a33  Slco3a1  Apc2  Slc7a2  Slc12a6  Slc6a17  Npc1  Slc23a2  Hiatl1  Slc31a1  Hiat1  Slc25a11  Slc44a2  Tusc3  Slc1a5  Slc25a3  Sfxn1 [1] "Sodium_channels"
[1] "TRP_channels"
5 genes:  Trpc3  Trpc4  Trpv2  Pkd2  Trpm7 

transporter_dictionary_for_table <- transporter_dictionary
for (i in names(transporter_dictionary)) {
  transporter_dictionary_for_table[[i]] <- cat(paste0(transporter_dictionary[[i]]))
}

transporter_dictionary_for_figure <- transporter_dictionary
#transporter_dictionary_for_figure$GPCR <- transporter_dictionary$GPCR[1:17]

ion_channel_genes <- c(transporter_dictionary_for_figure$Calcium_channels, transporter_dictionary_for_figure$Potassium_channels, transporter_dictionary_for_figure$Sodium_channels, transporter_dictionary_for_figure$TRP_channels)
ion_channel_genes <- unique(ion_channel_genes)
ion_channel_genes <- sort_gene_by_ubiquity(combined_cell, ion_channel_genes, "x")
transporter_dictionary_for_figure$ion_channel_genes <- ion_channel_genes

receptor_genes <- c(transporter_dictionary$Amine_receptors, transporter_dictionary$Cholinergic_receptors,transporter_dictionary$Glutamate_receptors, transporter_dictionary$Glycine_receptors)
receptor_genes <- unique(receptor_genes)
receptor_genes <- sort_gene_by_ubiquity(combined_cell, receptor_genes, "x")
transporter_dictionary_for_figure$receptor_genes <- receptor_genes

pump_genes <- c(transporter_dictionary$ABC_transporters, transporter_dictionary$ATPase)
pump_genes <- unique(pump_genes)
pump_genes <- sort_gene_by_ubiquity(combined_cell, pump_genes, "x")
transporter_dictionary_for_figure$pump_genes <- pump_genes

generate_dotplot(combined_cell, ion_channel_genes, paste("pc_Ion_Channels.csv"))

generate_dotplot(combined_cell, receptor_genes, paste("pc_Receptors.csv"))

generate_dotplot(combined_cell, transporter_dictionary_for_figure$GPCR, paste("pc_GPCR.csv"))

generate_dotplot(combined_cell, pump_genes, paste("pc_Pump.csv"))

NA
NA

# intersction: gene_symbols_slc & brain.markers.sandra
slc_targets.sandra <- intersect(tolower(gene_symbols_slc), tolower(brain.markers.sandra$Gene))
filtered.brain.markers.sandra <- brain.markers.sandra %>% filter(tolower(Gene) %in% slc_targets.sandra)
slc_targets.sandra <- filtered.brain.markers.sandra$Gene
write.csv(filtered.brain.markers.sandra, "filtered_Brain_EC_markers_fc1_padj005.csv", row.names = TRUE)

# dot plot of the expressions
dot_plot <- DotPlot(combined_cell, features = slc_targets.sandra[1:27], group.by = "organ") +
  labs(x = "Organ", y = "Gene") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  coord_flip()
ggsave("dot_plot_brain.pdf", plot = dot_plot, width = 7, height = 10)

sc analysis

my36colors <- c('#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3', '#5F3D69', '#C5DEBA', '#58A4C3', '#E4C755', '#F7F398','#585658',
         '#AA9A59', '#E63863','#E95C59', '#E59CC4', '#AB3282', '#23452F', '#BD956A', '#8C549C', '#476D87',
         '#9FA3A8', '#E0D4CA', '#E39A35', '#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B',
         '#712820', '#DCC1DD', '#CCE0F5', '#CCC9E6', '#625D9E', '#68A180', '#3A6963',
         '#968175')

DimPlot(combined_cell, reduction = "umap", group.by = "organ", label = TRUE, repel = TRUE,cols = my36colors)

deg analysis

library(tidyverse)
library(RColorBrewer)
library(scales)
library(reshape2)
library(tidyverse)
library(harmony)
library(readxl)
library(readr)

# brain slc degs
brain_slc_degs <- read_csv("~/Library/CloudStorage/GoogleDrive-schen601@usc.edu/My Drive/transportor/deg_SLC/Brain.csv")

plot_gene <- brain_slc_degs[,c("cluster","...1")]
colnames(plot_gene)[colnames(plot_gene) == "...1"] <- "gene"
plot_gene$gene <- sub("\\..*", "", plot_gene$gene)

#plot_gene <- plot_gene[1:20,]

pct_thres <- 0.2
fc_thres <- 0.2
pv <- 0.05

brain_slc_degs_filtered <- brain_slc_degs[(brain_slc_degs$p_val<pv & abs(brain_slc_degs$avg_log2FC) > fc_thres),]


brain_slc_degs_filtered <- brain_slc_degs_filtered[(((brain_slc_degs_filtered$avg_log2FC > fc_thres)&(brain_slc_degs_filtered$pct.1>pct_thres))|((brain_slc_degs_filtered$avg_log2FC < -fc_thres)&(brain_slc_degs_filtered$pct.2>pct_thres))),]


plot_gene <- brain_slc_degs_filtered[,c("cluster","...1")]
colnames(plot_gene)[colnames(plot_gene) == "...1"] <- "gene"
plot_gene$gene <- sub("\\..*", "", plot_gene$gene)

violin

# Load required libraries
library(ggplot2)
library(Seurat)
library(gridExtra)

# Define output PDF file
pdf("brain_deg_violin_plot.pdf", width = 6, height = 6)  # Adjust width & height as needed

# Set up plot counter
plot_list <- list()
plot_count <- 0
plots_per_page <- 4  # 2*2 grid per page

# Generate violin plots
candidates <- plot_gene$gene
candidates <- c(deg_v_plot1,deg_v_plot2,deg_v_plot3)

for (gene in candidates) {
  if (gene %in% rownames(combined_cell)) {  # Ensure the gene exists in the data
    p <- VlnPlot(combined_cell, features = gene, group.by = "organ", pt.size = 0.1) +
      ggtitle(gene) +
      theme(plot.title = element_text(hjust = 0.5, size = 8))

    plot_list[[length(plot_list) + 1]] <- p
    plot_count <- plot_count + 1

    # When reaching 9 plots, print the page
    if (plot_count %% plots_per_page == 0 || gene == tail(candidates, 1)) {
      do.call(grid.arrange, c(plot_list, ncol = 2, nrow = 2))
      plot_list <- list()  # Reset plot list
    }
  }
}

# Close the PDF file
dev.off()

break pdf

# Load necessary library
library(pdftools)

# Define input PDF file
input_pdf <- "raw_violin_plots.pdf"  # Replace with your actual PDF file

# Read the number of pages in the PDF
total_pages <- pdf_info(input_pdf)$pages  # Should be 26

# Define output file names
output_files <- c("raw_violin_part1.pdf", "raw_violin_part2.pdf", "raw_violin_part3.pdf", "raw_violin_part4.pdf")

# Define page splits (first 6, next 6, next 6, last 8)
page_splits <- list(1:6, 7:12, 13:18, 19:26)

# Loop through each subset and save to a new PDF
for (i in seq_along(output_files)) {
  pdf_subset(input_pdf, pages = page_splits[[i]], output = output_files[i])
  message("Saved: ", output_files[i])
}

# self
deg_v_plot1 <- c("Slco1c1", "Slc22a8", "Slc38a3", "Slc38a5", "Slc1a1", "Slc19a3", "Stra6","Slc6a13","Slc38a5")

#shared with Testis
deg_v_plot2 <- c("Mfsd2a", "Slc7a5", "Slc16a4", "Mfsd7c","Slc35f2", "Slc19a3", "Mfsd7c", "Slc35f2")

#w/ bladder
deg_v_plot3 <- c("Slc38a11", "Slc5a5")

which candidate fully show up

# Load required library
library(dplyr)

# Define the list of candidate genes (replace with your actual list)
candidates <- c(deg_v_plot1, deg_v_plot2, deg_v_plot3)  # Example candidate genes

# Define the required organ conditions
required_organs <- c("Brain", "Bladder", "Liver", "Heart", "Intestine",
                     "Colon", "Testis", "Muscle", "Lung", "Spleen", "Kidney")

# Filter for candidate genes only
all.markers.filtered.by.candidates <- all.markers %>%
  filter(gene %in% candidates)

# Count how many unique organs each gene appears in
gene_organ_count <- all.markers.filtered.by.candidates %>%
  group_by(gene) %>%
  summarize(organs_present = list(unique(cluster)), count = n())

# Check which genes appear in **all 10 organs**
genes_fully_present <- gene_organ_count %>%
  filter(all(required_organs %in% organs_present)) %>%
  pull(gene)  # Extract the gene names

# Print results
print("Genes found in all 10 organs:")
print(genes_fully_present)

dot plot

# Load required libraries
library(ggplot2)
#library(dplyr)

# Transform p-value
deg_table <- all.markers.filtered.by.candidates %>%
  mutate(p_size = -log10(p_val+0.0000000000000000000000000000000000000000000000000000000000000000000000000000000000000001))  # Transform P-value

# Create the dot plot
dot_plot <- ggplot(deg_table, aes(x = cluster, y = gene)) +
  geom_point(aes(size = p_size, color = avg_log2FC)) +  # Dot size = -log10(p_val), Color = log2FC
  scale_size_continuous(range = c(1, 8)) +  # Adjust dot size range
  scale_color_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0) +  # Color gradient
  theme_minimal() +
  labs(x = "Organ", y = "Gene", size = "-log10(P-value)", color = "Log2 Fold Change") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Display the plot
print(dot_plot)

Brain markers

heatmap

DoHeatmap(object = combined_hm, features = features, slot = "data", group.colors = my36colors)+ scale_fill_gradientn(colors = c('#3b69b5',  "white", '#b53b3b'))

#plot_gene <- read_excel("plotgene.xlsx")
plot_gene$cluster=factor(plot_gene$cluster,levels = unique(plot_gene$cluster))
plot_gene <- plot_gene%>%arrange(cluster,gene)

combined_cell@meta.data$cell_type=factor(combined_cell@meta.data$cell_type,levels = levels(plot_gene$cluster))

combined_cell@meta.data$CB=rownames(combined_cell@meta.data)
#integrated@meta.data=integrated@meta.data%>%inner_join(integrated,by="CB")
rownames(combined_cell@meta.data)=combined_cell@meta.data$CB

color_ct=c(brewer.pal(12, "Set3")[-c(2,3,9,12)],"#b3b3b3",
           brewer.pal(5, "Set1")[2],
           brewer.pal(3, "Dark2")[1],
           "#fc4e2a","#fb9a99","#f781bf","#e7298a")
names(color_ct)=levels(plot_gene$cluster)

### 主代码 ######################################################################
#vln.df=as.data.frame(integrated[["RNA"]]@data[plot_gene$gene,])
vln.df=as.data.frame(as.data.frame(as.matrix(GetAssayData(combined_cell[["RNA"]], slot = "data")))[plot_gene$gene, ])
vln.df$gene=rownames(vln.df)
vln.df=melt(vln.df,id="gene")
colnames(vln.df)[c(2,3)]=c("CB","exp")

anno=combined_cell@meta.data[,c("CB","cell_type")]
vln.df=inner_join(vln.df,anno,by="CB")
vln.df$gene=factor(vln.df$gene,levels = plot_gene$gene)
vln.df <- vln.df[!is.na(vln.df$cell_type), ]

# 当你想竖直方向排版这张图片时：
vln.df%>%ggplot(aes(cell_type,exp))+
  geom_violin(aes(fill=cell_type),scale = "width")+ 
  #如果想最终呈现出来的图，是根据基因涂色的，也就是一个基因一种颜色，应改为：aes(fill=gene)
  #一般而言，基因数多于细胞类型数，当根据基因涂色时，配色方案有点繁琐，所以不推荐aes(fill=gene)
  facet_grid(gene~.,scales = "free_y")+
  scale_fill_manual(values = color_ct)+
  scale_y_continuous(expand = c(0,0))+
  theme_bw()+
  theme(
    panel.grid = element_blank(),
    
    axis.title.x.bottom = element_blank(),
    axis.ticks.x.bottom = element_blank(),
    axis.text.x.bottom = element_text(angle = 45,hjust = 1,vjust = NULL,color = "black",size = 14),
    axis.title.y.left = element_blank(),
    axis.ticks.y.left = element_blank(),
    axis.text.y.left = element_blank(),
    
    legend.position = "none",
    
    panel.spacing.y = unit(0, "cm"),
    strip.text.y = element_text(angle=0,size = 14,hjust = 0),
    strip.background.y = element_blank()
  )
ggsave("brain_slc_markers.pdf",device = "pdf",width = 18,height = 90,units = "cm")